home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / listings / v_11_05 / 1105068a < prev    next >
Text File  |  1993-03-07  |  2KB  |  59 lines
  1. #include <math.h>
  2. #define SWAP(a,b) tempr=a;a=b;b=tempr
  3. void four1(float *data, int *nn, int *isign)
  4. {                /* altered for consistency
  5.                  * with original FORTRAN */
  6.     /* Press, Flannery, Teukolsky, Vettering "Numerical
  7.      * Recipes in C" tuned up ; Code works only when *nn is
  8.      * a power of 2  */
  9.     int n, mmax, m, j, i;
  10.     double wtemp, wr, wpr, wpi, wi, theta, wpin;
  11.     float tempr, tempi, datar, datai;
  12.     n = *nn * 2;
  13.     j = 0;
  14.     for (i = 0; i < n; i += 2) {
  15.     if (j > i) {        /* could use j>i+1 to help
  16.                  * compiler analysis */
  17.         SWAP(data[j], data[i]);
  18.         SWAP(data[j + 1], data[i + 1]);
  19.     }
  20.     m = *nn;
  21.     while (m >= 2 && j >= m) {
  22.         j -= m;
  23.         m >>= 1;
  24.     }
  25.     j += m;
  26.     }
  27.     theta = 3.141592653589795 * .5;
  28.     if (*isign < 0)
  29.     theta = -theta;
  30.     wpin = 0;            /* sin(+-PI) */
  31.     for (mmax = 2; n > mmax; mmax *= 2) {
  32.     wpi = wpin;
  33.     wpin = sin(theta);
  34.     wpr = 1 - wpin * wpin - wpin * wpin;    /* cos(theta*2) */
  35.     theta *= .5;
  36.     wr = 1;
  37.     wi = 0;
  38.     for (m = 0; m < mmax; m += 2) {
  39.         for (i = m; i < n; i += mmax * 2) {
  40.         j = i + mmax;
  41.         /* mixed precision not significantly more
  42.          * accurate here; if removing float casts,
  43.          * tempr and tempi should be double */
  44.         tempr = (float) wr *data[j] - (float) wi *data[j + 1];
  45.         tempi = (float) wr *data[j + 1] + (float) wi *data[j];
  46.         /* don't expect compiler to analyze j >
  47.          * i+1; original code stored data[j..]
  48.          * first and avoided using data[ri] */
  49.         data[i] = (datar = data[i]) + tempr;
  50.         data[i + 1] = (datai = data[i + 1]) + tempi;
  51.         data[j] = datar - tempr;
  52.         data[j + 1] = datai - tempi;
  53.         }
  54.         wr = (wtemp = wr) * wpr - wi * wpi;
  55.         wi = wtemp * wpi + wi * wpr;
  56.     }
  57.     }
  58. }
  59.